DIPGEOPR_03

class: inverse,left, middle
background-image: url(background.png)
background-size: cover

##Módulo 3: Percepción remota, análisis masivo y Google Earth Engine.

### Pre-procesamiento satelital I

Gabriel Castro

mail gabriel.castro.b@pucv.cl</a>

.large[<a href="https://www.pucv.cl/uuaa/site/edic/base/port/labgrs.html">LabGRS</a> | Octubre 2024]

---
class: center,middle
background-image: url(data:image/png;base64,#labgrs_logo.png)
background-size: 35%

---
## Contenidos

.pull-left[
- Imagenes comerciales: Planet Scope Skysat.

- Combinaciones de bandas.

- Calibración radiométrica. (TOA radianza - TOA reflectancia)

- Corrección atmosférica.

- Filtros de calidad de la información satelital.

- Compensación topográfica por iluminación.
 
 ]
 
.pull-right[
<img src="data:image/png;base64,#https://raw.githubusercontent.com/allisonhorst/stats-illustrations/main/rstats-artwork/r_rollercoaster.png" width="650px"/>
© Allison Horst
]

---
### Recordemos.

---
### 1) Imagenes satelitales PlanetScope - Skysat.

#### - [**Información general**](https://assets.planet.com/docs/Planet_Combined_Imagery_Product_Specs_letter_screen.pdf)
#### - Especificaciones tecnicas. 
#### - Niveles de procesamiento. 
#### - Metadatos e interpretación de la información.

---

---
<center>
<img src="data:image/png;base64,#planet_catalogo.PNG" alt="drawing" width="950" height="500"/>
<center>

[**Especificaciones catálogo Skysat**](https://developers.planet.com/docs/data/skysat/)

---
### 2) Combinaciones de bandas.

A partir de las signaturas espectrales es posible entender el comportamiento de los objetos frente a diferentes longitudes de onda.

Esto nos permite crear combinaciones de bandas en los canales **RGB** para resaltar elementos específicos en superficie.

```r
#librerias necesarias ---- 
library(satellite)
library(tidyverse)
library(terra)
library(raster)

#Producto Analytic - ortho. 
skysat <- rast("Skysat_ViñadelMar/SkySat_08_feb_2024_vina.tif")

par(mfrow = c(1, 3))
plotRGB(skysat, r=3, g=2, b=1, stretch = "lin", axes = TRUE)
plotRGB(skysat, r=4, g=2, b=1, stretch = "lin", axes = TRUE)
plotRGB(skysat, r=3, g=4, b=1, stretch = "lin", axes = TRUE)
dev.off()
```

---
### 2) Combinaciones de bandas (Skysat).

<center>
<img src="data:image/png;base64,#comb_bandas_skysat.PNG" alt="drawing" width="1000" height="500"/>
<center>
---

### 2) Combinaciones de bandas (Landsat).

---
### 3) Calibración radiométrica.

---
### 3) Calibración radiométrica.

En un gran número de sensores (TM, ETM+, OLI, etc.) se establece la conversión a **radianza espectral** en base a una **relación lineal** existente entre el nivel digital bruto y los valores de radianza (unidades físicas). Sin embargo, **siempre es importante** revisar los manuales de usuario de cada producto con el que estemos trabajando.

De tal forma que la ecuación que describe la radianza espectral, según Chander et al. (2009) es:

---
### 3) Calibración radiométrica (skysat)

---
### 3) Calibración radiométrica (skysat)

---
### 3) Calibración radiométrica (skysat) 
#### Niveles digitales a radianza (TOA).

```r
# 5.- Calibración radiometrica ND TOA radianza (manual) ---- 
TOA_rad_skysat <- skysat * 0.01 # factor escalado

# Comparar histogramas:
par(mfrow = c(1, 2))
hist(skysat[[1]], main = "histograma banda azul ND")
hist(TOA_rad_skysat[[1]], main = "histograma banda azul TOA")
dev.off()
```

---

<center>
<img src="data:image/png;base64,#ND_TOA_REF.PNG" alt="drawing" width="1050" height="600"/>
<center>
---
### 3) Calibración radiométrica 
#### Niveles digitales a reflectancia aparente (TOA reflectancia)

- Transformar la radianza a reflectancia al tope de la atmosfera es un paso intermedio a la corrección atmosférica final.

- Si bien la mayoría de los métodos de corrección transforman automáticamente los niveles digitales a reflectancia superficial, es importante conocer el **porcentaje** de energía que es reflejada por los objetos en la superficie y que finalmente es captada por el sensor.

- Al ser la reflectancia al tope de la atmosfera, se le conoce como **reflectancia aparente**.
---
### 3) Calibración radiométrica 
#### Niveles digitales a reflectancia aparente (TOA reflectancia)

---
### 3) Calibración radiométrica 
#### Niveles digitales a reflectancia aparente (TOA reflectancia).

```r
# 6.- Calibración radiometrica ND TOA reflectancia ----
## proceso manual --------------------------------------
d <- calcEarthSunDist(date = "2024-02-08", formula = "ESA") # AU
ESUN <- c(2009.28, 1820.25, 1583.3, 1114.22) # (W * m2 * sr1 * um1)
pi <- pi # 3.14
## transformar angulos a radianes ----------------------
sun_elev <- 65.4*pi/180 # Sun elev radianes 
sun_azimuth <- 312.9*pi/180 # Sun azimuth radianes 
sun_zenith <- (90 - 65.4)*pi/180 # solar zenith

TOA_reflectance_skysat <- (pi*TOA_rad_skysat*d^2)/(ESUN*cos(90-sun_elev))
```

---

---

### 4) Corrección atmosferica.

Clasificación de correcciones atmosféricas:

- **Mediciones in situ**

- **Mixtas** Ejemplo: Landsat Ecosystem Disturbance Adaptative Processing System (LEPADS) y Land Surface Reflectance Code (LaSRC)

- **A partir de modelos físicos de transferencia radiativa** Ejemplos: Second Simulation of a Satellite Signal in the Solar Spectrum (6S vector version), Multi-angle implementation of atmospheri correction (MAIAC).

- **Con datos de la propia imagen** Ejemplos: Simple Dark object Substraction (sDOS), Sen2cor, ACOLITE.
---
<center>
<img src="data:image/png;base64,#reF_app_SR.PNG" alt="drawing" width="1000" height="600"/>
<center>

---
### 4) Corrección atmosferica.
<center>
<img src="data:image/png;base64,#irradiacion.PNG" alt="drawing" width="700" height="500"/>
<center>
---
### 4) Corrección atmosferica.
<center>
<img src="data:image/png;base64,#softwares.PNG" alt="drawing" width="900" height="500"/>
<center>

---
### 4) Corrección atmosferica. 
 Métodos darkest object substraction:

En la literatura podemos encontrar diferentes métodos DOS los cuales se basan en principios similares.

El paradigma inicial señala que algunos pixeles son cuerpos oscuros cuya reflectancia seria cercana a 0. Esto indica que la radianza captada por el sensor de estos pixeles oscuros, estaría influida por el efecto de la **bruma**.

Para estos métodos se estima que los piexeles oscuros reflejan al menos entre el 1-2% de la energía. (0.01)

El método **darkest pixel subtraction** busca reducir o eliminar este efecto de la bruma para obtener la reflectancia superficial.
---
### 4) Corrección atmosferica.
<center>
<img src="data:image/png;base64,#solarAngles.PNG" alt="drawing" width="920" height="510"/>
<center>

---
### 4) Corrección atmosferica. 
<center>
<img src="data:image/png;base64,#DIFF_METODOS2.PNG" alt="drawing" width="960" height="560"/>
<center>

---
<center>
<img src="data:image/png;base64,#METODOS.PNG" alt="drawing" width="730" height="650"/>
<center>
---
### Metodo DOS2 Chavez 1996.

Este método de corrección por DOS presento una mejora respecto al **simple darkest object subtraction** al tener en consideración la corrección de la transmitancia atmosférica desde la superficie terrestre hacia el sol.

También conocido como método CosTz, esta mejora utiliza el ángulo zenital del sol el cual está asociado a la transmitancia atmosférica para las fechas especificas en las cuales se tomó la imagen.

Sin embargo este método sigue considerando la Transmitancia de la energía del sol a la tierra (Tv) como constante (no pierde energía en el viaje de entrada se mantiene con valor 1 como el caso del sDOS).

---
### Comparemos ambos metodos: sDOS 1989 v/s DOS2 1996. 
#### Parametros iniciales sDOS:

```r
## Corrección atmosferica metodo sDOS ---------- 
Tz <- 1 
Tv <- 1
Edown <- 0 
## ----------------------------------------------
skysat_raster <- skysat %>% stack()
B <- calcDODN(skysat_raster[[1]]) 
G <- calcDODN(skysat_raster[[2]])
R <- calcDODN(skysat_raster[[3]])
NIR <- calcDODN(skysat_raster[[4]])
DNm <- c(B,G,R,NIR)
hazeValues <- DNm*0.01 
```
**Funcion calcDODN del paquete satellite calcula de forma automatica los valores de radiancia minimo para cada banda considerando n cantidad de piexeles (20% aproximadamente)**

---
### Calculo parametros finales sDOS 1989:

##### 1) L_001: Calculo de la reflectancia minima del 1%.
##### 2) SHV_rad: valores de bruma para cada banda considerando la influencia atmosferica.
##### 3) L_haze: calculo de la influencia atmosferica en la radianza del sensor.
##### 4) SR_sDOS: Calculo sDOS.

```r
L_001 <- 0.01 * (ESUN[1:4] * cos(sun_zenith) * Tz + Edown) * Tv / (pi * d^2)
SHV_rad <- hazeValues

L_haze <- SHV_rad[1:4] - L_001[1:4]

SR_sDOS <- pi * (TOA_rad_skysat[[1:4]] - L_haze[1:4]) * d^2 / 
 (Tv * ESUN[1:4] * cos(sun_zenith) * TAUz) + Edown

plotRGB(SR_sDOS, r=3, g=2, b=1, stretch = "lin", axes = TRUE, 
        main = "Reflectancia superficial metodo sDOS")

writeRaster(SR_sDOS, filename = "Reflectancia_superficial_metodo_sDOS.tif", overwrite = T)
```
---
### Calculo parametros finales DOS2 1996:

#### 1) Los parametros son los mismos donde lo unico que cambia es Tz.

```r
## Corrección atmosferica metodo sDOS ---------- 
Tz <- cos(sun_zenith)
Tv <- 1
Edown <- 0 
## ----------------------------------------------

L_001 <- 0.01 * (ESUN[1:4] * cos(sun_zenith) * Tz + Edown) * Tv / (pi * d^2)
SHV_rad <- hazeValues

L_haze <- SHV_rad[1:4] - L_001[1:4]

SR_DOS2 <- pi * (TOA_rad_skysat[[1:4]] - L_haze[1:4]) * d^2 / 
 (Tv * ESUN[1:4] * cos(sun_zenith) * TAUz) + Edown

plotRGB(SR_DOS2, r=3, g=2, b=1, stretch = "lin", axes = TRUE, 
        main = "Rflectancia superficial metodo DOS2")

writeRaster(SR_DOS2, filename = "Reflectancia_superficial_metodo_DOS2.tif", overwrite = T)
```
---
### Comparativa.

```r
### AED resultados obtenidos
par(mfrow = c(2,3))
boxplot(TOA_reflectance_skysat[[2]], outline = F, main = "Boxplot TOA_ref", 
        col = "green", ylim = c(0, 0.30))

boxplot(SR_sDOS[[2]], outline = F, main = "Boxplot sDOS", 
        col = "lightblue", ylim = c(0, 0.30))

boxplot(SR_DOS2[[2]], outline = F, main = "Boxplot DOS2", 
        col = "red", ylim = c(0, 0.30))
# ----------------------------------------------------
hist(TOA_reflectance_skysat[[2]])
hist(SR_sDOS[[2]])
hist(SR_DOS2[[2]])
dev.off()
```

---
### Comparativa.
<center>
<img src="data:image/png;base64,#comparacion_metodos.PNG" alt="drawing" width="1000" height="550"/>
<center>
---
### Resumén calibración radiometrica y corrección atmosférica.

- Ya sea estemos trabajando con colecciones satelitales Landsat, Sentinel, Planet, Worldview etc. **Siempre se deben revisar los manuales de uso de cada producto**. Si al momento de recibir las imagenes estas no cuentan con un manual de uso, es obligatorio exigirlos.

- **Siempre revisar los metadatos** de los productos. Es importante tener en cuenta cuales son los parámetros que tenemos y cuales deben ser estimados por otros medios. Esto permitirá definir la línea de procesamiento a seguir.

- Siempre evaluar que los resultados obtenidos tengan **coherencia**. Valores sobre o bajo el rango esperado, problemas en la visualización, pueden obedecer a errores en la línea de procesamiento. **Es importante realizar un pequeño AED** para garantizar la integridad de los datos y eliminar artefactos.

---
### Bibliografía.

Thomas Nauss, F.D. (2021) satellite: Handling and Manipulating Remote Sensing Data, Satellite - classes and methods for satellite data.

Chavez, Jr, Pat. (1996). Image-Based Atmospheric Corrections - Revisited and Improved. Photogrammetric Engineering and Remote Sensing. 62. 1025-1036.

Chavez & P.S, 2015. Image-Based Atmospheric Corrections - Revisited and Improved. , no. January 1996.

Moran, M.S., Jackson, R.D., Slater, P.N. y Teiuet, P.M., 1992. Evaluation of Simplified Procedures for Retrieval of Land Surface Reflectance Factors from Satellite Sensor Output. , vol. 184, pp. 169-184.

Song, C., Woodcock, C.E., Seto, K.C., Lenney, M.P. y Macomber, S.A., 2000. Classification and Change Detection Using Landsat TM Data : When and How to Correct Atmospheric Effects ? , vol. 4257, no. 00.

Wang, Y. , Wang, X. y He, H., 2019. An Improved Dark Object Subtraction Method for Atmospheric Correction of Remote Sensing Images [en línea]. S.l.: Springer Singapore. ISBN 9789811399176. Disponible en: http://dx.doi.org/10.1007/978-981-13-9917-6_41.

---
class: inverse middle
![logo](data:image/png;base64,#LABGRS-PUCV.png)